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Abstract 

Currently,  the  746th  Test  Squadrons  (746th  TS)  Central  Inertial  and  GPS  Test 
Facility  (CIGTF)  operates  one  of  the  most  accurate  truth  reference  systems,  called 
the  CIGTF  Reference  System  (CRS).  CIGTF  will  be  replacing  the  CRS  with  a  new 
references  system  called  UHARS  (Ultra  High  Accuracy  Reference  System).  UHARS 
will  differ  from  CRS  by  adding  the  ability  to  use  a  non-GPS  pseudolite  system,  as  a 
new  measurement  source.  This  research  effort  describes  the  design  of  the  extended 
Kalman  filter  which  is  developed  in  AFIT’s  SPIDER  filter  framework  which  imple¬ 
ments  a  tightly-coupled  pseudolite/INS  filter. 
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INITIAL  IMPLEMENTATION  AND  TESTING  OF  A 
TIGHTLY-COUPLED  IMU/PSEUDOLITE  SYSTEM 

I.  Introduction 


1.1  Motivation 

The  Global  Positioning  System  (GPS)  is  widely  accepted  as  a  valuable  tool  for 
position,  navigation,  and  timing.  As  the  system  continues  to  proliferate,  the  threats 
to  the  system  continue  to  likewise  grow.  This  creates  a  need  to  develop  alternative 
navigation  systems  that  can  continue  to  operate  in  a  GPS  degraded  or  denied  envi¬ 
ronments.  This  generates  the  requirement  to  develop  a  system  that  operates  for  test 
systems  in  GPS  degraded  environments.  One  such  alternative  to  relying  on  GPS  for 
a  truth  reference  is  the  currently  emerging  LUtra  High  Accuracy  Reference  Systems 
(LIHARS)  [11].  The  UHARS  project  is  currently  under  development  by  the  746th 
Test  Squadron  Central  Inertial  and  GPS  Test  Facility  (CIGTF)  and  will  be  further 
discussed  in  Section  1.2.  In  the  past,  to  test  GPS  under  jammed  conditions,  only  a 
single  frequency  was  jammed  in  order  to  leave  the  remaining  frequency  to  obtain  a 
truth  solution  [28].  At  other  times,  both  frequencies  would  be  jammed,  and  then  the 
truth  relics  on  a  smoothed  solution  based  on  the  inertial  navigation  system  (INS). 
This  system  includes  the  use  of  ground-based  transceivers,  or  pseudolites,  to  trans¬ 
mit  the  signals  used  to  develop  the  navigation  solution.  By  incorporating  pseudolites, 
which  operate  on  a  different  radio- frequency  (RF)  band,  more  realistic  test  scenarios 
can  occur  where  all  GPS  RF  bands  are  jammed  or  degraded. 
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1.2  Background 


The  746th  Test  Squadron  Central  Inertial  and  GPS  Test  Facility  (CIGTF)  part¬ 
nered  with  Locata  and  AFIT  to  develop  an  improved  truth  reference  system.  Shockley 
provides  a  linage  of  the  navigation  solution  truth  reference  systems  used  by  the  746th 
TS  and  Figure  1.1  illustrates  the  reference  systems  development  [23].  Previously, 
CIGTF  employed  other  reference  systems  to  assess  position,  velocity,  and  time.  The 
first  system  the  CIGTF  developed  was  the  Completely  Integrated  Reference  Instru¬ 
mentation  System  (CIRIS)  [18].  CIRIS  was  a  Fortran  based  system,  that  depended 
on  an  INS  with  beacon  measurements.  The  CIGTF  replaced  CIRIS  in  1990  with 
the  CIGTF  High  Accuracy  Post-processing  Reference  System  (CHAPS)  [26].  This 
reference  system  was  developed  in  1990  and  began  to  include  GPS  signals [23].  The 
CIGTF  Reference  System  (CRS),  replaced  CHAPS  in  and  is  the  current  system  [18]. 
CIGTF  is  near  the  completion  of  replacing  the  CIGTF  Reference  System  (CRS)  with 
the  Ultra  High  Accuracy  Reference  System  (UHARS)  currently  under  development 
as  the  next  generation  truth  reference  system. 


1960's 
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k _ _ _ . _ 


C - \ 
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Figure  1.1.  746th  Test  Squadron  Truth  Reference  Systems 
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1.3  Problem  Definition 


This  research  seeks  to  implement  an  algorithm  to  process  pseudolite  measurements 
tightly-coupled  with  an  inertial  measurement  unit  (IMU)  or  an  inertial  navigation  sys¬ 
tem  (INS)  for  obtaining  a  post-processed  navigation  solution.  The  pseudolite  solution 
is  expected  to  obtain  a  similar  level  of  accuracy  of  GPS  carrier-phase  solutions [11]. 
The  development  of  this  new  truth  reference  system  must  ultimately  have  the  accu¬ 
racy  of  at  least  the  current  systems  accuracy  under  dynamic  conditions.  The  UHARS 
set-up,  further  discussed  in  Chapter  II  and  Chapter  III,  is  capable  for  measuring  po¬ 
sitions  for  both  ground  and  flight  testing[23].  Data  collected  in  early-November  2014 
flight  tests  are  used  for  this  research.  For  the  developed  system,  carrier-phase  based 
navigation  solutions  were  generated  using  the  Sensor  Processing  and  Inertial  Dynam¬ 
ics  Error  Reduction  (SPIDER)  software  written  in  M  AT  LAB® .  This  framework, 
discussed  in  Chapter  III,  allows  for  an  expedited  implementation  of  a  Kalman  filter. 

1.4  Proposed  Solution 

In  this  research,  measurements  from  the  CRS  and  UHARS  systems  are  processed 
and  the  results  are  analyzed.  A  tightly-coupled  Locata  Loclite/IMU  navigation  solu¬ 
tion  with  cycle-slip  detection  and  a  state  reset  algorithm  is  developed  and  compared 
against  the  truth  from  the  Novatel  truth  reference  software  Waypoint.  The  solution  is 
analyzed  against  the  truth  to  determine  the  accuracy  of  the  developed  algorithm.  To 
implement  the  algorithm,  the  MAT  LAB®- based  SPIDER  framework  is  employed. 
The  results  of  the  algorithm  are  further  discussed  in  Chapter  IV. 

SPIDER. 

The  Sensor  Processing  and  Inertial  Dynamics  Error  Reduction  (SPIDER)  pack¬ 
age  for  M  AT  LAB® ,  developed  by  AFIT,  allows  for  quick  implementation  of  new 
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navigation  filter  algorithms  by  having  an  established  modular  framework  in  place. 
Currently  the  SPIDER  package  allows  the  user  to  select  sensors,  and  add  scenario 
specific  inputs  as  needed  to  a  Kalman  filter  [25].  During  this  research  effort  a  GPS  dif¬ 
ferential  carrier-phase  solution  and  a  Locata  pseudolite  carrier-phase  model  (including 
a  pseudolite  tropospheric  delay  model)  arc  implemented  in  SPIDER  and  tested  with 
inputs  from  test  flight  data. 

1.5  Thesis  Overview 

Background  information  and  a  survey  of  literature  is  provided  in  Chapter  II.  The 
algorithms  to  model  the  pseudolite  navigation  solution  and  implemented  tropospheric 
models  are  discuses  in  Chapter  III.  The  processed  results  for  the  flight  and  ground 
tests  are  analyzed  in  Chapter  IV.  The  results  and  discussion  of  the  likely  sources  of 
error  are  analyzed  as  well.  Chapter  V  suggests  future  work  needed  and  summarizes 
relevance  of  this  research. 
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II.  Background 


2.1  Introduction 

This  chapter  presents  the  background  information  for  the  development  of  the  nav¬ 
igation  filter,  failure  detection/rejection  algorithm,  and  the  troposphere  delay  models 
for  this  research.  In  conducting  the  general  outline  of  the  research  approach,  previous 
works  from  AFIT  were  used  [3]  [7]  [10]. 

2.2  Related  Works 

In  2004,  Bouska  developed  simulations  that  studied  the  possibility  that  pseudo- 
lites  can  measure  position  accuracy  to  a  0.10  to  0.15  meter  level  [7].  Carrier-phase 
measurements  were  used  in  the  simulation  using  floating-point  ambiguities.  Bouska’s 
simulations  focused  on  a  pseudolite-only  solution  and  were  modeled  using  first-order 
Gauss-Markov  functions  to  simulate  both  noise  and  uncertainty.  Bouska  stated  in  his 
research  that  centimeter-level  accuracy  is  possible  with  a  pseudolite  only  reference 
system.  Further,  a  tropospheric  scale  factor  error,  which  would  appear  to  look  like 
constant  biases  in  the  pseudolite  measurements.  This  errors  can  be  filtered  by  adding 
a  state  or  states  that  would  scale  the  tropospheric  corrections  by  so  linear  factor  [8] . 
His  research  has  since  provided  a  cornerstone  for  other  research  on  pseudolite  based 
navigation  to  include  this  document. 

The  optimal  geometry  of  the  Locata  pseudolites  has  been  evaluated  by  several  re¬ 
searchers  [12]  [13]  [32],  Crawford  focused  on  using  the  geometric  normalized  accuracy 
of  precision  (GNAP)  rather  that  the  geometric  dilution  of  precision  (GDOP)  as  a 
metric  [12],  Crawford  showed  through  his  research  that  an  orbiting  pseudolite,  or 
an  inverted  pseudolite,  holds  promise  to  increase  the  level  of  precision.  Having  an 
orbiting  pseudolite  may  not  be  always  practical,  so  alternatively  he  showed  that  as 
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the  elevation  angle  increased,  the  geometric  errors  decreased.  Further,  accuracy  was 
also  improved  by  the  number  of  deployed  pseudolites[12],  Dai  further  discussed  the 
optimal  geometry  of  the  pseudolite  positions,  and  also  explored  an  inverted  pseudolite 
configuration  [13].  J.J.  Wang  discussed  the  near-far  problem  in  consideration  of  the 
optimal  design[32].  The  near  far  problem  is  depicted  in  Figure  2.2. 


The  near-far  problem,  shown  in  Figure  2.1  illustrates  the  useable  bounds  of  the 
pseudolites.  For  the  purpose  of  this  research  effort,  the  aircraft’s  flight  over  the  Loclite 
array  is  not  affected  by  the  near-far  problem.  In  2006  Amt  using  a  Locata  pseudolite 
only  reference  system  demonstrated  that  ground  and  air  vehicle  positions  could  be 
accurately  measured  [3] .  Amt  used  the  carrier-phase  measurements  and  pseudoranges 
to  estimate  the  position  solution.  The  pseudorange  measurements  applied  a  batch 
least  squares  algorithm,  while  the  carrier-phase  measurements  were  integrated  using 
floating  point  ambiguity  resolution  in  a  Kalman  filter.  Further,  the  development  of 
a  height  (altitude)  determination  for  the  pseudolite  measurements  contributed  to  the 
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overall  progress  of  implementing  the  UHARS  program.  Shockley  presented  simula¬ 
tion  and  real  data  results  for  a  Locata  pseudolite  positioning  system  that  used  single 
differenced  carrier-phase  measurements  [3].  His  research  demonstrated  that  a  pseudo¬ 
lite  positioning  system  is  capable  of  decimeter-level  accuracy  when  the  measurement 
errors  are  mitigated  through  modcling[4].  The  tropospheric  delay  error,  which  Shock- 
ley  chose  to  model  by  using  a  scale-factor  in  a  manner  similar  to  Bouska,  allowed  for 
results  with  improved  accuracy [24],  In  2009,  Ciampa  furthered  the  existing  research 
by  proposing  a  failure  detection  algorithm  for  Locata  pseudolites[10].  Ciampa  com¬ 
pared  failure  detection  residual  monitoring  algorithms  based  on  thresholding  and  a 
moving  window  approach.  After  completion  of  his  research,  the  recommendation 
was  to  use  a  combined  approach  of  residual  monitoring  within  the  framework  of  a 
Kalman  filter.  Further  discussion  of  residual  monitoring  is  discussed  in  Section  2.8. 
Several  researchers  have  performed  a  comparison  of  various  pseudolite  tropospheric 
delay  models  producing  similar  results[16]  [31].  In  Wang’s  research,  the  Radio  Tech¬ 
nical  Commission  for  Aeronautics  model  (RTCA),  modified  RTCA,  Bouska- Raquet, 
Saastamoinen,  Nicll  mapping  function  and  UNB3m  models  were  compared  by  sim¬ 
ulation  [31].  The  RTCA,  Hopheld  and  Saastamoinen  models  were  compared  by  So 
using  flight  test  data  with  radiosonde  meteorological  data  [16].  In  both  cases  no  single 
model’s  performance  had  the  best  overall  results. 

2.3  Navigation  Reference  Frames 

A  position  and  velocity  trajectories  are  related  in  relative  terms  by  using  a  ref¬ 
erence  frame.  Below  are  descriptions  of  reference  frame  coordinate  systems  that  are 
employed  in  this  research.  Figure  2.2  depicts  the  several  reference  frames  defined  for 
this  research[27]. 
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Earth  Centered  Earth  Fixed. 


The  ECEF  frame  has  the  origin  at  the  center  of  the  Earth.  The  z-axis  pints 
through  the  North  Pole,  the  x-axis  points  to  the  intersection  of  the  Equator  and 
Prime-Meridian  and  the  y- axis  points  orthogonal  to  both  the  x  and  z  axes  to  complete 
the  right-handed  coordinate  system  [27].  The  ECEF  frame  is  used  in  this  research 
when  it  is  desirable  to  have  all  axis  in  units  of  meters.  This  is  of  particular  interest 
when  finding  the  distance  between  to  objects.  Figure  2.2  provides  a  visual  of  the 
ECEF  reference  frame  along  with  the  navigation  reference  frame  [30]. 


Local 

Meridian 


Figure  2.2.  ECEF  and  Navigation  Frames.  Used  with  original  author’s  permission[30]. 


Inertial  Frame. 


The  o  inertial  frame  like  the  ECEF  has  an  origin  at  the  center  of  the  Earth.  The 
z-axis  of  the  inertial  frame  also  points  through  the  North  Pole.  The  x-axis,  points 
towards  first  star  of  Aries.  The  y-axis  like  in  the  ECEF  points  orthogonal  to  both 
the  x  and  z  axes  completing  the  right-handed  coordinate  system  [27]. 

Platform  or  Body  Frame. 

The  platform  or  body  frame  reference  frame  uses  the  aircraft  itself.  The  origin 
of  the  body  frame  located  at  the  center  of  mass  or  some  other  fixed  point  of  the 
aircraft.  The  x-axis  points  towards  the  nose  of  the  aircraft,  the  y-axis  points  towards 
the  right  wing  of  the  aircraft,  and  the  z-axis  point  orthogonal  to  the  x  and  y  axis 
downwards [2 7].  The  sensor  frame  is  the  frame  in  which  the  sensor  measurements  are 
directly  measured.  It  can  be  related  to  the  body  frame  through  the  use  of  lever  arms 
(for  translation)  and  also  a  relative  orientation  if  needed  [27]. 

Geodetic  Coordinate  system. 

The  geodetic  coordinate  system  is  a  widely  used  reference  frame.  The  axes  this 
coordinate  systems  uses  are  fixed  by  the  latitude  and  longitude  lines  of  the  Earth  and 
an  altitude  off  of  the  WGS-84  model[27]. 

2.4  Global  Positioning  System  (GPS)  Signals 

The  Global  Positioning  Systems  is  the  current  benchmark  of  high  precision  nav¬ 
igation  reference  systems.  GPS  satellites  are  identified  by  a  unique  pseudorandom 
noise  (PRN)  with  a  provided  a  navigation  message.  Within  the  navigation  message, 
the  satellites  constellation’s  ephemeris  is  given  so  that  the  position  of  the  satellite 
is  known  to  some  degree  of  error.  The  GPS  satellites  currently  can  transmit  3  fre- 
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quencies,  however,  two  of  the  frequencies,  often  call  the  legacy  frequencies,  are  used 
most  widely  for  navigation  solutions.  These  frequencies  are  the  LI  (which  operates  at 
1575.42  MHz)  and  L2  (which  operates  at  1227.6  MHz).  The  LI  frequency  has  both 
the  course  acquisition  (C/A)  and  precise  (P)  codes,  where  the  L2  frequency  only 
carriers  the  P  code.  A  GPS  position  solution  most  commonly  achieves  a  position 
solution  based  on  some  combination  of  the  following  three  measurements:  pseudor¬ 
ange,  carrier  phase,  and  Doppler.  The  pseudorange  and  carrier-phase  measurements 
are  further  discussed  below. 


Pseudorange  Measurements. 


The  pseudorange  measurements,  or  code  measurements,  are  based  on  the  differ¬ 
ence  between  the  time  of  the  signal  transmission  and  reception.  The  pseudorange 
measurements,  although  are  widely  available,  have  many  errors  associated  with  the 
measurements.  Some  of  these  errors  include  the  troposphere  delay  error,  ionosphere 
error,  clock  errors,  measurement  noise,  multipath  and  variety  of  other  errors  factors 
[15]  [17].  Since  several  errors  exist,  the  possibility  exists  that  even  a  well  modeled 
measurement  will  achieve  a  navigation  solution  for  an  air  vehicle  with  approximately 
upwards  of  10  meters  of  error  [14].  The  measurements  for  pseudorange  is  expressed  as 
2.1. 


p  =  r  +  c(Stm  -  5ts )  +  T  +  mp  +  vp 


(2.1) 


where: 

p  =  GPS  pseudorange  measurement  ( meters ) 
r  =  Range  from  the  user’s  receiver  to  the  transmitter  (meters) 
c  =  speed  of  light  in  a  vacuum  (meters / second) 

8tm  =  User  clock  error  (seconds) 

8tp  =  Transmitter  clock  error  (seconds) 
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T  =  pseudorange  error  due  to  Tropospheric  delay 
I  =  pseudorange  error  due  to  Ionospheric  delay  ( meters ) 
rrip  =  pseudorange  error  due  to  multipath  ( meters ) 
vp  =  pseudorange  error  clue  to  receiver  noise  ( meters ) 

Carrier-Phase  Measurements. 

The  carrier-phase  measurements  are  more  accurate  than  pseudorange  measure¬ 
ments  since  the  noise  and  multipath  are  much  smaller.  Carrier-phase  measurements 
are  based  on  the  number  of  cycles  of  the  carrier  frequency  over  the  distance  between 
the  satellite  and  receiver.  An  integer  ambiguity  state,  N,  is  created  and  becomes  a 
source  of  error  for  the  carrier-phase  measurements  until  it  is  resolved.  Once  the  inte¬ 
ger  ambiguity,  is  resolved,  the  integer  ambiguity  may  suffer  from  cycle  slips.  When  a 
cycle  slip  occurs  the  integer  ambiguity  typically  is  reset  and  must  again  reinitialize. 
The  equation  for  carrier-phase  measurements  is  expressed  below  in  Equation  2.2 [IT] . 

(j)  =  A-1(r  +  c(8tm  -  8tp )  +  T  +  +  v^)  +  N  (2.2) 


where: 

0  =  carrier-phase  measurement  ( cycles ) 

A  =  carrier-phase  wavelength  (cycles /second) 
nip,  =  carrier-phase  error  due  to  multipath  ( meters ) 

Vtj,  =  carrier-phase  error  due  to  receiver  noise  (meters) 
N  =  carrier-phase  integer  ambiguity  ( cycles ) 
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2.5  Pseudolites  for  Navigation 


Interest  in  augmenting  position,  navigation  and  timing  (PNT)  infrastructures  with 
pseudolites  has  been  increasing  [6].  When  compared  against  GNSS,  pseudolite  are 
cost-efficient,  relatively  easy  to  held,  and  can  provide  the  same  level  or  better  posi¬ 
tioning  accuracy  [11].  Further,  pseudolites  can  operate  on  separate  radio  frequency 
bands  than  GNSS  and  therefore  can  operate  in  GNSS  jammed  environments  [5]  [6]. 
The  pseudolites  used  in  this  research  are  Locata  Loclites.  The  partnership  between 
the  746th  Test  Squadron  at  Holloman  AFB  and  the  Locata  Corporation  has  provided 
an  opportunity  to  preform  flight  testing  by  the  746th  TS  culminating  into  the  next 
generation  flight  reference  system.  Research  from  AFIT  has  made  significant  contri¬ 
butions  to  this  framework  since  2004  [21]  [3]  [10]  [24],  Figure  2.3  illustrates  that  the 
pseudolites  operate  on  a  ’master-slave’  system.  The  master  pseudolite  time  synchro¬ 
nizes  all  the  slave  pseudolite  for  precision  timing. 
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Figure  2.3.  Locata  Loclie  Master-Slave  Configuration 

The  pseudolites  systems  can  generate  the  same  observables  as  GPS;  however, 
they  are  often  ground-based  systems.  The  pseudolites  pseudorange  and  carrier-phase 
measurements  at  preformed  just  as  the  GPS  measurement  as  explained  above  in 
Equation  2.1  and  Equation  2.2.  One  significant  change  is  that  in  the  place  of  the 
GPS  pseudorandom  noise  (PRN)  as  a  method  to  identify  the  space  vehicles,  a  naming 
convention  was  created  for  the  pseudolites.  The  term,  LRPN,  is  used  for  ease  and 
commonality  between  GPS  and  the  pseudolites. 
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2.6  Kalman  Filter 


The  Kalman  filter  is  a  linear,  optimal,  recursive  estimator  that  provides  estimate 
over  time.  The  estimates  can  either  exist  in  continuous  time  or  discrete  time;  how¬ 
ever,  for  computational  ease,  the  Kalman  filter  usually  exist  in  a  discrete  fixed-time 
interval.  The  use  of  Kalman  filtering  is  widely  used  for  error-state  navigation  filters, 
particularly  the  extended  Kalman  filter.  This  section  provides  a  brief  description  of 
the  discrete  Kalman  filter  algorithm [19]. 

The  Discrete  Kalman  Filter  Model. 

For  aviation,  the  Kalman  filter  often  begins  a  known  location  of  the  aircraft. 
As  the  aircraft  begins  to  move,  the  dynamics  of  the  aircraft  are  propagated  until  a 
measurement  from  the  INS/IMU  or  positioning  system  updates  the  location.  The 
Kalman  filter  has  two  main  components:  the  time  update  phase,  or  propagate,  and 
the  measurement  correction  phase,  or  update  [9].  Additive  white  noise  is  introduced  in 
both  the  propagation  phase  and  during  the  measurement  update.  The  propagation 
phase  of  the  Kalman  filter  occurs  at  some  defined  interval  in  which  the  modeled 
systems  dynamics  are  projected  forward  to  the  time  of  the  next  measurement.  The 
state  estimate  at  this  time  provides  the  a  priori  estimate  for  the  next  measurement 
incorporation.  When  the  next  measurement  is  available,  the  filter  The  Figure  2.4 
depicts  the  cycle  of  the  Kalman  filter. 
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Time  Update  +  Process  Noise 
(Propagation) 


Figure  2.4.  Kalamn  Filter  Cycle 

By  applying  the  linear  differential  equation  characterized  in  linear  time,  is  shown 
equation  2.3 [19]  .  The  input  noise  of  the  system’s  covariance  is  shown  in  equation  2.4 
[19] 

x(t)  =  Fx(t)  +  Gw(t)  (2.3) 

Q (t)  =  E  w(t)w(t)T  (2.4) 

The  system’s  dynamics  that  are  modeled  with  states  are  represented  by  the  Fx(t)  in 
Equation  2.3.  The  system  noise  is  represented  by  the  Gw(t)  term,  where  G  is  the 
system  noise  mappi 

Pxx(t)  =  E  x(t)x(t)T  (2.5) 

Pxx(t)  =  FP  xx(t)  +  Pxx(t)FT  +  GQ(t)GT  (2.6) 

Where  P  represents  the  state  covariance  matrix  and  the  remaining  terms  have  been 
previously  defined.  The  next  set  of  equations  that  the  discrete  KF  usess  are  the 
equations  to  bring  the  filter  into  a  discrete  form.  The  variable  /c,  represents  the 
current  discrete  time.  In  Equation  2.7  ,  the  state  transition  matrix  is  brought  into 
a  discrete  form  represented  by  the  symbol  $[19]  .  The  noise  also  in  brought  into  a 
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discrete  form  and  is  represented  by  the  Qd  represented  by  Equation  2.8. 


$  =  eF(At)k 


(2.7) 


Qd 


E 


(2.8) 


The  next  steps  in  the  discrete  KF  is  to  apply  the  discrete  state  propagation  and  the 
discrete  state  covariance  propagation  shown  in  equations  2.9  and  2.10  . 


xfc  =  <£xfc_i  +  wfc  (2.9) 

Pfc  =  $Pfc_!$T  +  Qd  (2.10) 

The  discrete  state  estimate,  Xfc  and  the  discrete  covariance  matrix,  variance  Pfc  is 
propagated  forward  in  time.  The  previous  estimated  state  mean  and  covariance  ma¬ 
trix,  represented  by  xa,_i  and  P^-i  respectively,  are  applied  to  the  Equations  2.11 
and  2.12  .  These  equations  then  calculate  the  estimated  state  vector,  and  covariance 
represented  by  x^  and  P]7  . 

+  (2.n) 

P*  =  #Pt1*T  +  Qd  (2.12) 

Once  a  measurement  becomes  available  to  the  Kalman  filter,  the  measurement  up¬ 
dates  can  take  place.  The  z*,  term  shows  the  filter’s  estimation  of  the  desired  state [19]. 
The  H^Xfc  term  represents  where  the  state  estimate  is  transformed  into  measurement 
space  by  the  sensitivity  matrix,  H/.,  of  a  sensor,  and  the  relationship  to  the  states[19]. 
The  measurement  noise,  v^,  is  typically  defined  by  white  Gaussian  noise.  Lastly  the 
R  fc,  or  the  measurement  noise  covariance  matrix  or  measurement  uncertainty  covari¬ 
ance  matrix,  is  defined  by  the  type  of  sensor/measurement  being  used [9].  Equations 
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2.14  and  2.14  show  how  the  linear  measurement  model  f  the  KF  and  the  noise  covari¬ 


ance  matrix  are  calculated. 


R/ 


E 


T 


(2.13) 

(2.14) 


The  update  phase  of  the  KF  is  the  last  part  of  the  KF  algorithm.  This  phase,  the 
measurement  estimates,  z/,.  ,  are  inputted  into  the  algorithm  to  bring  about  a  linear 
optimal  estimate  of  the  state  estimates  and  state  covariance  estimates.  The  KF 
estimate  is  founded  from  a  minimum  mean  square  estimation  error [9].  Equation  2.15 
shows  the  Bayesian  optimal  estimate  of 


Xjt  =  E(xfc|zfc) 


(2.15) 


The  sensors  measurements  are  feed  into  the  Kalman  filter  at  some  time  interval  defined 
by  each  sensor.  At  this  time  interval  a  measure  update  occurs  [9].  The  Kalman  gain, 
equation  show  in  2.16,  a  posteriori  estimate  covariance,  equation  shown  in  2.17,  and 
a  posteriori  state  estimate  are  calculated,  equation  shown  in  2.18. 


14  =  P1THl’(H,P^Hf  +  R)-1  (2.16) 

p*  =  (I  -  A'<:Ht)Pt:  (2.17) 

X;7  -  !!*'>.(.  - H,.xk  /  (2.18) 

2.7  Extended  Kalman  Filters 

The  extended  Kalman  filter  (EKF)  is  the  extension  of  the  Kalman  filter  that  is  able 
to  estimate  nonlinear  state  dynamics  or  measurements.  The  extended  Kalman  Filter 


17 


(EKF)  is  similar  to  the  KF,  in  that  the  EKF  also  a  recursive  estimator  that  provide 
time  based  updates  that  is  linearized  around  a  estimate [9]  [20].  The  extended  Kalman 
filter  unlike  the  Kalman  filter,  is  sub-optimal,  except  when  estimating  linear  problems, 
in  which  case  is  identical  to  the  Kalman  filter [20] .  The  EKF  approaches  nonlinearities 
by  linearizing  the  state  dynamics  or  measurement  around  a  point.  Each  iteration  the 
is  linearized  by  a  first-order  Taylor  series  approximation  of  the  nonlinearity.  Further 
discussion  of  the  EKF  is  found  in  [20]. 

Extended  Kalman  Filters  Equations. 

The  propagation  phase  for  the  EKF  consists  of  similar  equations  to  that  of  the 
KF.  The  EKF  propagate  phase  is  much  like  that  of  of  the  KF;  however,  the  state 
estimate,  f(xfc),  can  be  nonlinear.  During  the  update  phase  updated  by  a  non-linear 
function  represented  by  f(xfc).  The  residual,  sensitivity  matrix,  and  Kalman  gain 
differs  from  the  KF  shown  in  Equations  2.19,  and  2.20. 

+KtK  “  h(*i))  (2-19) 

P t  =  (I  -  KtHt)Pt-  (2.20) 

The  EKF  measurement  model,  shown  in  Equation  2.21.  is  designed  to  handle  nonlin¬ 
earities.  The  state  observabilities  are  mapped  by  the  sensitivity  matrix  Hx*,  within 
the  KF  algorithm  and  is  replaced  by  passing  the  estimated  states  through  a  non¬ 
linear  function  h(xfc)[9][20].  The  non-linear  measurement  function  h(xfc)  that  passes 
through  the  EKF,  is  approximately  linearized  within  the  measurement  matrix  by  the 
use  of  a  first-order  shown  in  Equation  2.22  Taylor  series  expansion. 


z  k  =  h(xfc)  +  vfc 


(2.21) 
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(2.22) 


d 

Hfc  =  7— h(xfc) 

<9xfc 

2.8  Filtering  Techniques 

When  running  a  EKF  in  MATLAB,  processing  time  and  numerical  instability  may 
become  an  issue.  Techniques  developed  to  counter  some  of  these  obstacles  that  arise 
is  by  using  the  upper  diagonal  matrix  or  the  inverse  covariance  matrix [20]. 

Error  States  Model. 

Rather  that  modeling  the  states  of  the  actual  measurement  of  the  system’s  dy¬ 
namics,  it  is  advantageous  to  model  the  error  states.  Error  state  models  are  common 
techniques  to  model  stochastic  process,  and  are  widely  used  for  navigation.  The 
error  state  are  the  simply  put  the  difference  between  the  ’truth’  and  the  estimate. 
From  error  state  models  the  truth,  estimate,  and  error  can  be  easily  compared  for 
analysis  [9]. 

Residual  Monitoring. 

The  integrity  of  a  Kalman  filter  is  entirely  dependent  on  the  ability  to  correctly 
model  the  system.  One  metric  of  a  navigation  filter’s  performance  is  the  magnitude 
of  the  residuals.  This  thesis  the  residuals  are  monitored  for  both  the  carrier  phase 
and  pseudorange  measurements  and  discussed  in  Chapter  IV [20]. 

Upper  Diagonal  Factorization  Matrix. 

The  U-D  factorization  implementation  of  the  Kalman  filter  involves  the  use  of  an 
upper  triangular  matrix  U  and  a  diagonal  matrix  D  defind  in  Equation  2.23 [9] . 

P  =  UDUt  (2.23) 
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Equation  2.23  is  called  the  U-D  factorization  of  P.  This  allows  the  propagation  and 
update  of  U  and  D  instead  of  P.  decreasing  the  numerical  issues  in  the  implemen¬ 
tation  of  the  Kalman  filter’s  algorithm.  The  variation  of  the  matrix  elements  gets 
compressed  so  that  the  dynamic  range  is  reduce.  This  results  is  higher  efficiency  to 
process  the  algorithm.  The  initial  conditions  start  with  x(fd)  and  P(f,+),  and  then 
finds  U(t+)  and  D (fd)  using  U-D  factorization.  The  propagation  operation  for  the 
states  follows  the  same  as  the  KF. 

Inverse  Covariance  Matrix. 

The  inverse  covariance  matrix  is  used  when  the  values  of  the  covariance  become 
large.  This  technique  is  often  used  in  optimal  smoothing  discussed  below  in  Section 
2.8. 


Optimal  Smoothers. 

Optimal  smoothing  consists  of  three  sub-types:  fixed-point  smoothing,  fixed- 
interval  smoothing,  and  fixed-lag  smoothing[20].  The  algorithm  for  fixed-point  smooth¬ 
ing  runs  with  selected  points  backwards  in  time.  In  fixed-interval  smoothing,  the 
algorithm  run  with  a  fixed  time  interval  using  a  posteriori  measurements.  Fixed-lag 
smoothing,  would  be  used  if  a  near-real  time  filter  is  desired  and  the  algorithm  run 
on  a  delayed  fixed  time  interval  so  that  the  real-time  system  can  still  estimate  in  real¬ 
time.  In  this  research  fixed-interval  smoothing  is  the  most  appropriate  smoothing 
algorithm  for  the  post  processed  data  analysis[20].  A  fixed- interval  Kalman  smoother 
is  divided  into  two  parts,  with  the  first  part  consisting  of  a  standard  KF  running 
forward  in  time  from  the  initial  time  to  to  the  final  time  t/  and  the  second  part  with 
a  filter  running  backward  in  time  from  tf  to  to-  Upon  completion  of  the  forward 
and  backward  filter  runs,  the  state  estimates  are  optimally  combined  to  create  a 
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“smoother”  estimate.  For  the  forward  filter,  a  priori  information  about  x(to)  is  used 
to  initialize  the  KF,  while  for  the  backward  filter,  the  “initial”  condition  is  established 
with  x(tf)  as  a  random  vector  with  no  available  a  priori  statistical  information.  A 
fixed-interval  KF  smoother  requires  one  to  estimate  the  states  both  forward  and 
backward  in  time[20].  When  implementing  a  backward  KF,  the  covariance  of  the 
state  at  the  final  time  (first  time  index  for  the  backward  method)  is  initialized  to 
be  infinite,  meaning  there  is  no  information  available  about  the  states  of  the  system. 
The  derivation  of  the  backward  KF,  using  the  inverse  covariance  form  discussed  in 
Section  2.8  and  the  optimal  smoothing  of  the  two  KFs  is  found  below.  It  begins  by 
initializing  the  backward  state  estimate  and  the  inverse  covariance  matrix  to  zero [20]. 
The  fixed-interval  Kalman  smoother  will  produce  the  most  optimal  result  possible[9]. 
This  method  is  only  valid  in  a  post-test  analysis  since  it  requires  all  measurements 
from  time  initial  to  time  final.  This  method  requires  more  computational  resources 
than  the  standard  KF  implementation,  but  in  the  post-test  analysis  environment  this 
increase  in  computational  overhead  is  not  a  significant  concern.  Kalman  filters  are 
recursive  in  nature,  and  run  in  forward  time.  The  information  that  is  taken  into 
account  currently  in  a  Kalman  filter  is  always  from  the  previous  time  epoch. 


Navigation 

Solution 


Figure  2.5.  Reference  Frames 
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2.9  Troposphere  Delay  Models 


For  navigation  employing  the  use  of  pseudolites,  the  tropospheric  delay  becomes 
a  major  source  of  errors[ll].  To  achieve  sub-decimeter  level  accuracy,  the  tropo¬ 
spheric  model  must  be  an  appropriate  estimation  of  the  error  induced  upon  the  signal. 
The  troposphere  is  a  persistent  source  of  error  when  propagating  a  radio  frequency 
through  this  medium.  Since  the  troposphere  is  a  neutral,  non-dispersive  medium, 
the  error  caused  by  the  delay  of  the  signal  can  be  modeled.  To  appropriately  model 
the  error,  the  troposphere  must  further  be  broken  down  into  two  components:  the 
hydrostatic  component,  more  commonly  referred  to  the  dry  component;  and  the  non¬ 
hydrostatic  component,  which  is  more  commonly  referred  to  as  the  wet  component  of 
the  troposphere[31].  The  dry  component  of  the  tropospheric  error  model  consists  of 
80-90%  of  the  error.  The  dry  component  is  predictable  based  on  the  season  (position 
of  the  Earth),  altitude  and  general  location  of  the  receiver.  Thus  the  dry  component 
can  be  modeled  in  such  a  way  to  ensure  that  99%  of  the  error  can  be  modeled  for 
correctly.  On  the  other  hand,  the  wet  component  of  the  tropospheric  error  model  only 
consists  of  10-20%  of  the  error,  and  generally  is  difficult  to  predict.  It  is  from  the 
difficultly  of  modeling  the  wet  component  of  the  troposphere  that  is  the  significant 
differentiator  that  separates  the  various  troposphere  delay  models [16].  Several  tropo¬ 
spheric  delay  models  that  are  suitable  for  pseudolite  positioning  have  been  developed, 
however,  the  models  fall  into  three  categories.  The  first  category  is  the  integration 
method.  An  example  of  this  method  is  the  Hopheld  model.  The  second  category  is 
this  differencing  method  such  as  the  modified  RTCA  model.  The  third  category  is  the 
scaling  method.  In  the  scaling  method,  a  scaling  factor  is  added  to  the  entire  model 
in  the  state  dynamics  matrix.  A  comparison  of  these  methods  will  be  compared  in 
this  research  and  will  be  further  discussed  in  Chapter  III. 
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2.10  Summary 


The  information  provided  in  this  chapter  provides  the  general  overview  of  the 
signals  and  methods  that  will  be  used  in  development  of  the  research.  GPS  and 
pseudolite  signals,  both  having  pseudorange  and  carrier  measurements  are  inputs 
into  the  extended  Kalman  filter  as  will  be  presented  in  Chapter  III. 
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III.  Navigation  Algorithm 


3.1  Overview 

This  chapter  lays  out  the  development  of  the  algorithms  that  implement  the 
tightly-coupled  pseudolite/INS  integration.  Descriptions  of  the  state-models,  mea¬ 
surement  models,  tropospheric  model,  and  other  necessary  equations  are  contained 
in  this  chapter.  The  navigation  filter  is  first  examined. 

Navigation  Filter. 

The  navigation  filter  implemented  in  this  research  is  derived  from  an  extended 
Kalman  filter,  discussed  in  Chapter  2,  implemented  in  the  SPIDER  framework,  dis¬ 
cussed  in  Chapter  1.  The  filter  employs  an  error-state  approach  while  tightly  coupling 
the  pseudolite  measurements  to  the  inertial  measurements  to  reduce  error.  Although, 
the  pseudolites  measurement  data  is  received  at  a  10  Hz  rate,  the  measurements  can 
be  down  sampled  to  any  rate  desired. 

Inertial  Model. 

As  discussed  in  Chapter  II,  state  space  modeling  is  a  common  method  used  for 
Kalman  filtering.  An  error  state  model  is  developed  to  include  all  states  of  interest. 
By  modeling  the  error  states,  rather  than  estimating  the  actual  position,  velocity,  tilt 
angle,  clock,  and  so  on,  this  filter  will  rather  output  the  errors  of  these  states,  this 
can  then  be  applied  to  the  inertial  solution  to  obtain  the  final  estimate. 

State  Dynamics  Model. 

The  model  used  in  this  research  is  the  Pinson-15.  As  stated  in  Chapter  II  [30], 
this  model  is  an  error  states  model  of  an  IMU  or  INS.  This  research  employes  previ- 
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ously  used  model,  the  states  consist  of  a  three  axis  position,  velocity,  tilt,  gyro  bias, 
and  acceleration  bias [30].  These  are  computed  in  the  navigation  reference  frame,  as 
defined  in  Chapter  II  2.3.  The  locations  of  the  individual  states  are  listed  below  to 
form  the  first  15-states  of  the  state  dynamics  model,  or  the  F  matrix. 

X\  :  Latitude  (North)  Position  Error  (rad) 
x2  :  Longitude(East) Position  Error  (rad) 

X3  :  Altitude  (Down)  Position  Error  (m) 

X\  :  North  Velocity  Error  (m/s) 
x5  :  East  Velocity  Error  (m/s) 

Xq  :  Down  Velocity  Error  (m/s) 

X7  :  North  Tilt  Angle  Error  (rad) 
x$  :  East  Tilt  Angle  Error  (rad) 

Xg  :  Down  Tilt  Angle  Error  (rad) 

X\q  :  X-axis  Accel  Bias/Error  (m/s2) 
x\\  :  Y-axis  Accel  Bias/Error  (m/s2) 

X12  ■  Z-axis  Accel  Bias/Error  (m/s2) 

Xi3  :  X-axis  Gyro  Bias/Error  (rad/ sec) 
xu  :  Y-axis  Gyro  Bias/Error  (rad/ sec) 

X15  :  Z-axis  Gyro  Bias/Error  (rad/ sec) 

The  state  dynamic  model,  Fins  will  form  initially  as  a  15x15  matrix  as  shown  below 
in  Equation  3.1.  As  other  states  are  added,  including  clock  error,  and  ambiguities, 
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the  state  dynamics  model  will  grow  while  maintaining  an  n  —  by  —  n  matrix [30]. 
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(3.1) 


where 

C® :  navigation  to  Earth  frame  rotation  matrix 
C^:body  to  navigation  frame  rotation  matrix 
Cg :  Earth  to  navigation  frame  rotation  matrix 
fn:  force  in  the  navigation  frame 
G:  noise  intensity  matrix 
Vte\  Earth  rate 

Ta:  time  constant  for  accelerometer  bias 
Tu:  time  constant  for  accelerometer  bias 


In  addition  to  the  states  in  the  dynamics  matrix  is  an  accompanying  system  noise 
intensity  matrix,  or  the  G  matrix.  The  noise  intensity  matrix  describes  the  growth 
in  the  uncertainty  over  time  and  is  dependent  on  the  IMU/INS  used  for  propagation. 
In  this  research  an  enhanced  embedded  global  positioning  system  inertial  navigation 
system (EEGI)  is  used.  From  the  G  matrix  values,  each  states  value  will  vary  as  the 
system  propagates  forward.  Ideally,  the  G  matrixs  values  correctly  model  the  EEGIs 
actual  noise  characteristics.  The  Equations  3.2  and  3.3  removed  the  biases  in  the  gyro 
and  accelerometer  noise,  while  Equation  3.4  models  the  drift  in  the  accelerometer  and 
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gyro  [30], 


ffns  —  fbMS  +  abias  +  w{NS 

Uibins  =  Ubias  +  U ib  +  WINS 

^ bias  &bias  “1“  bias 

where 

fins\  forces  on  the  force  measurements  in  the  INS 
fbias :  bias  in  the  forces 
afejas:  acceleration  bias 

vJjns '■  white  guassian  noise  acting  upon  the  forces 
uj\b ,,, :  gyro  forces  in  the  body  frame 
Ubias'-  gyro  bias 

w°ins '■  white  guassian  noise  on  the  gyro  measnrments 
wbias'-  white  guassian  noiise  of  the  acceleration  bias 


(3.2) 

(3.3) 

(3.4) 


The  resulting  initial  noise  intensity  matrix  results  in  the  follow  shown  in  Equation 
3.5.  With  the  following  DCMs  defined  in  literature  and  identity  matrices  [27]. 
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Clock  Error. 


Clock  error  states  are  important  when  using  both  GPS  and  pseudolite  measure¬ 
ments  in  order  to  reduce  error.  Two-states  are  added  to  the  F  and  Q  matrixes  for  the 
clock  errors,  one  representing  the  clock  bias,  while  the  other  represents  the  derivative 
of  the  clock  bias  or  the  clock  drift.  These  states  will  be  used  directly  in  the  measure¬ 
ment  models  of  the  GPS  or  pseudolite  models  [15].  The  clock  errors  are  modeled  in 
the  filter  by  Equation  3.6 [29]. 

&clk\  0  1  %clk  i  ^  'cl  k  i 

=  +  (3.6) 

&clk2  0  0  *£c/fc2  UJrlk2 

GPS  and  Pseudolites. 

The  GPS  error-states  added  will  depend  on  which  measurement  method  is  being 
used.  For  pseudoranges,  the  clock  error  state  model  suffices.  Carrier  phase  measure¬ 
ments  do  require  additional  states  to  estimate  the  ambiguities.  Each  carrier-phase 
estimate  requires  an  additional  state  to  estimate  the  cycle  ambiguities.  This  states 
also  provides  insight  into  cycle  slip  detection  and  reset  for  updating  the  measurements 
in  the  case  that  a  cycle  slip  occurs  [26]  [29]. 

Final  State  Dynamics  And  Noise  Intensity. 

Now  that  as  the  applicable  states  are  created,  these  matrices  are  used  to  propagate 
the  filter  forward  in  time  by  using  the  equations  laid  out  previously  in  Chapter  II. 
The  dynamics  matrix  varies  with  each  epoch  depending  on  the  number  of  phase 
ambiguities  tracked,  but  will  consist  of  at  least  17  rows  and  columns  consisting  of 
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Pinson  15  terms,  and  the  clock  terms[30]. 


7  IN  Sib 
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(3.7) 


The  augmented  state  vector  that  is  derived  from  the  dynamics  matrix  consists  the 
diagonal  of  that  matrix.  The  noise  matrix  consists  of  an  n  —  by  —  m  matrix,  where 
n  is  equal  to  the  number  of  states  in  the  F  matrix,  and  m  is  equal  to  the  number  of 
measurements  included. 


3.2  Measurement  Model 

The  measurement  model  is  derived  using  the  measurements  available  to  the  filter 
at  a  moment  in  time.  In  this  research,  the  pseudoranges  and  carrier-phase  measure¬ 
ments  are  commonly  used.  The  h,  and  linearized  H.  matrices  contains  all  values 
derived  from  the  measurements  that  will  update  the  filter  with  the  equations  as  pre¬ 
viously  introduced  in  Chapter  2. 

Ultimately,  the  H  matrix  will  consist  of  an  m  —  by  —  n  matrix  where  n  is  the 
number  of  states  and  m  is  equal  to  the  number  of  measurements.  In  this  filter  the 
measurement  matrix [30].  The  measurement  uncertainty  or  measurement  noise  are 
modeled  in  the  R  matrix.  These  values,  like  the  values  of  the  G  matrix  provides 
the  appropriate  noise  when  correctly  modeled.  All  the  errors  are  applied  and  the 
corrected  measurements  are  non-linear.  Since  these  measurements  are  non-linear, 
for  each  epoch  the  measurement  model  is  linearized  using  a  first  order  Taylor  series 
approximation  [30].  The  pseudoranges  uses  the  estimates  from  the  clock  error  states 
and  the  location  of  the  pseudolites  to  provide  measurement  updates  to  the  filters 
position  estimate  of  the  mobile  receiver.  The  final  equation  to  derive  a  measurement 
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from  a  pseudorange  input  is  shown  in  Equation  3.8 [30]. 


Pn 


zu)2}  +  cSt 


(3.8) 


The  pseudorange  postion  estimate,  p,  is  base  on  the  xn,  yn,  and  zn  transmitter  loca¬ 
tions  differenced  from  the  estimated  mobile  receiver  location  ( xu ,  yu,  and  zu).  The 
cSt  term  represents  the  filter’s  estimated  clock  error  in  meters.  The  carrier  phase 
measurements  will  carry  addition  states  to  estimate  the  ambiguities,  described  by  the 
variable  Nn  in  units  of  cycles.  Depending  on  the  frequency  of  the  systems  used,  the 
cycles  are  converted  into  meters,  by  dividing  the  number  of  cycles  by  the  wavelength 
of  the  frequency  [29]. 


T- \j  Uxi,  -  x„)2  +  (y„  -  i y2  +  (z  -  iuf]  +  Nn+c5t 


(3.9) 


The  carrier-phase  position  estimate,  <f>n,  differs  from  the  pseudorange  position  esti¬ 
mate  by  including  a  term  to  convert  to  units  of  cycles  by  dividing  by  the  wavelength 
of  the  carrier,  An,  and  adding  in  the  ambiguity  term,  Nn.  The  resulting  h  vector 
will  appear  with  a  combination  of  the  pseudorange  measurements  and  carrier-phase 
measurements  from  Equations  3. 10 [29]. 


h [x(t  )  = 


{(xn  -  x uf  +  (yn  -  yu)2  +  (z  -  zu )2}  +  cSt 


t~  y  {(xn  —  x. u)  +  (yn  -  yu)  +  (z  -  zu)  }  +  Nn  +  c5t 


(3.10) 


The  H  matrix  will  appear  with  a  combination  of  the  pseudorange  measurements 
and  carrier-phase  measurements  after  taking  the  Jacobian  derivatives  from  Equations 
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shown  in  Equation  3. 11  [29]. 
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This  matrix  was  then  converted  to  the  latitude,  longitude,  and  altitude  states 
using  the  approach  given  in  [22],  Note  that  the  l’s  in  the  last  columns  of  the  H 
matrix  correspond  to  the  ambiguity  states  associated  with  each  measurement. 


3.3  Lever  Arms 

When  the  navigation  solution  is  calculated,  it  is  often  generated  at  the  location 
of  the  sensor.  In  this  research,  the  navigation  solutions  are  generated  at  the  inertial 
unit.  The  pointing  vector  from  a  sensor  to  the  inertial  unit  is  required  to  shift  the 
navigation  solution  to  a  common  point.  This  is  accomplish  by  using  direction  cosine 
matrices  and  pointing  vectors[30].  Figure  3.1  depicts  the  lever  arms  from  an  IMU  to 
an  RF  receiver. 
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Figure  3.1.  Example  of  Lever  Arms  from  an  Antennas  to  the  Interial 

3.4  SPIDER  Filter  Sequence 

The  EKF  will  run  through  the  SPIDER  framework  in  a  multitude  of  steps  [25]. 
First  the  specific  parameters  are  inputted  into  the  software.  The  parameters  include, 
but  are  not  limited  to,  the  timeframe  to  run  the  filter,  the  sensors  selected  for  the 
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measurement  updates,  noise  values,  lever  arms,  and  tuning  values.  The  states  are 
then  initialized  following  the  model  described  in  Section  3.2.  The  remaining  filter 
states  needed  for  the  sensors  are  then  initialized.  These  filter  states  include  the  clock 
terms,  barometer,  and  carrier-phase  ambiguities.  The  measurements  then  update 
the  Kalman  filer  and  and  update  the  navigation  solution.  The  output  is  then  saved, 
and  will  continue  for  the  remaining  measurement  epochs.  If  needed  a  final  lever  arm 
correction  is  applied  to  shift  the  solution  to  a  desired  location. 


Figure  3.2.  Spider  Filter  Sequence 
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3.5  Pseudolite  Tropospheric  Delay  Models 

Tropospheric  delay  as  discussed  in  related  works  [11]  [31],  is  a  dominant  error 
source  for  pseudolite  based  measurements.  When  calculating  the  tropospheric  delay, 
the  refractivity  for  the  hydrostatic  (dry)  and  non-hydrostatic  (wet)  components  must 
first  be  calculated  [31]. 

Ndry  =  77.6-^  (3.12) 

Where: 

Nr)ry  is  the  dry  refraction  index  of  the  medium 
P  is  the  pressure  at  the  pseudolite  measured  in  mmHG 
T  is  the  temperature  in  degrees  Kelvin 

7.4475(T  -  273 

Nwet  =  22770  •  -4  •  10  T  —  38.3  (3-13) 

1  - 

Where: 

Nwet  is  the  wet  refraction  index 
/  is  the  relative  humidity  at  the  pseudolite 
T  is  the  temperature  in  degrees  Kelvin 

Although  there  exist  several  tropospheric  models  for  tropospheric  delay,  three 
models  were  selected.  The  selected  models  were:  a  ray-tracing  modified  Hopheld 
model,  the  Bouska-Rauet  model  and  an  scaling  model.  These  models  at  most  only 
requires  inputs  of  the  altitude  of  the  pseudolite,  the  model  receiver,  elevation  angle, 
and  the  distance  between  the  pseudolite  and  mobile  receiver [16]  [8]. 

The  modified  Hopheld  troposphere  delay  model  is  applied  in  a  ray-tracing  format 
is  applied  by  calculating  the  altitude  segment  that  the  signal  passes  through.  At  each 
altitude  segment,  a  delay  term  is  estimated,  and  the  sum  of  these  terms  are  the  final 
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estimate  of  the  total  tropospheric  delay  terms  for  that  measurement.  An  illustration 
of  this  is  seen  in  3.3.  The  Hopheld  model  developed  in  1969  has  long  been  regarded 
as  a  well  performing  tropospheric  delay  model.  The  Hopheld  model  despite  being  one 
of  the  earliest  differencing  models  is  widely  used  still  for  GNSS.  The  equation  for  the 
Hopheld  model  is  given  in  3.14  and  3.15[31] . 


Ntrop(h)  =  NlJ 


K  -  h 
h* 
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(3.14) 


and 


Kop  =  K r 


h*  —  h 

K 


10  6ds 


(3.15) 


Where: 

Nlrop:iefT  activity 
h:  height 

Wf^irefr  activity  at  surface 
h*:scale  height 
A(rop:tropospheric  delay 
s:  signal  transmission  path 

Ray  tracing  allows  the  ability  to  segment  the  tropospheric  delay  area  into  smaller 
areas  in  order  to  obtain  a  better  result  [16].  This  is  done  by  diving  the  distance  the 
height,  h,  of  the  mobile  receiver  to  the  pseudolite  into  n  many  segments.  These 
segments  will  also  correspond  to  a  distance,  rn,  between  the  transmitter  and  the 
mobile  receiver.  Then  the  rn  and  hn  terms  feed  into  the  tropospheric  delay  models 
replacing  the  distance,  R,  and  height,  h,  terms. 
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Figure  3.3.  Ray  tracing  for  Tropospheric  Delay 


The  Bouska  and  Raquet  model  was  developed  based  on  the  RTCA  model  from 
an  internal  memo  for  Van  Dierendonck.  This  differencing  method  tropospheric  delay 
model  was  developed  in  2003  specifically  for  pseudolites.  The  equation  is  found  in 
3.16. 

Atrop  =  2  •  10-5AhRu- 
Where: 

Atrop  Tropospheric  delay 
N (surface  refractivity 

Ru:  range  from  pseudolite  to  reference  receiver 
h*)0:scale  height 

hs: height  of  the  reference  receiver 
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Ray  tracing  allows  the  ability  to  segment  the  tropospheric  delay  area  into  smaller 
areas  in  order  to  obtain  a  better  result.  This  is  done  by  diving  the  distance  the 
height,  H ,  of  the  mobile  receiver  to  the  pseudolite  into  n  many  segments.  These 
segments  will  also  correspond  to  a  distance,  rn,  between  the  transmitter  and  the 
mobile  receiver.  Then  the  rn  and  hn  terms  feed  into  the  tropospheric  delay  models 
replacing  the  distance,  R,  and  height,  H,  terms. 

The  index  of  refraction  model  is  the  scaling  method  of  calculating  the  tropospheric 
delay  [1] .  The  method  apples  a  scale  factor,  and  is  given  below  in  Equations  3.17,  3.18, 
and  3.19. 


Atrop  =  ND  ■  10“6 

(3.17) 

Nwet  =  71.2962^  +  375463  ^ 

(3.18) 

Ndry  =  77.689  (P°^eo) 

(3.19) 

3.6  Carrier  Phase  Cycle  Slip  Detection  and  Repair 

A  cycle  slip  detection  and  repair  method  is  important  for  when  correcting  carrier 
phase  measurements  [3].  Multipath  can  cause  a  poor  measurement  resulting  in  a  cycle 
slip.  The  other  common  cause  of  cycle  slips  in  this  data  is  that  a  phase  loop  lock 
is  lost  during  the  measurement  epoch  and  thus  causes  the  loss  of  that  measurement, 
resulting  in  a  cycle  slip  as  defines  in  the  developed  algorithm.  The  cycle-slip  detection 
algorithm  is  implemented  using  the  current  and  previous  epochs  measured  data  [29] . 
By  differencing  the  current  and  previous  carrier  phase  measurement  and  then  applying 
a  mean  Doppler  correction  term  the  slip  threshold  term  is  calculated.  A  carrier-phase 
threshold  is  set,  and  if  the  delta  carrier-phase  is  found  to  be  above  the  slip  threshold, 
a  cycle  slip  was  declared,  further  referred  to  as  a  threshold  slip.  Further,  cycle  slips 
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are  declare  if  two  consecutive  measurements  were  not  present,  further  referred  to  as 
a  time  slip. 

Slip  =  <t>t-  4>t- 1  +  [( Dopt  +  Dopt- 1)/2]  •  St  (3.20) 

Where 

c f>t  =  carrier-phase  measurement  at  current  epoch  (cycles) 

4>t-i  =  carrier-phase  measurement  at  previous  epoch  (cycles / second) 

Dopt  =  Doppler  measurement  at  current  epoch(meters) 

Dopt-i  =  Doppler  measurement  at  previous  epoch(mefers) 

St  =  Time  in  seconds  between  epochs  (meters) 

In  the  case  of  threshold  and  time  slips,  the  reinitialization  process  takes  place.  This 
process  consists  of  reinitializing  the  carrier-phase  ambiguity  state.  The  reinitialize 
state  in  this  algorithm  is  calculated  by  differencing  the  carrier  phase  measurement 
with  the  pseudorange  measurement  converted  to  cycles.  Further,  the  off-diagonal 
covariance  terms  of  the  pseudolite  are  reset  to  zero [15].  The  new  ambiguity  estimate 
is  recalculated  based  of  the  difference  between  the  carrier  phase  measurement  and  the 
estimated  receiver  position. 

In  the  case  of  threshold  and  time  slips,  the  repair  process  takes  place.  The  repair 
process  consists  of  re-initializing  the  carrier-phase  ambiguity  state.  The  reinitialize 
state  in  this  algorithm  is  calculated  by  differencing  the  carrier  phase  measurement 
with  the  filter  estimated  position  converted  to  cycles.  Further,  the  covariance  terms 
are  reset  to  zero  for  that  ambiguity  and  the  variance  term  is  also  reinitialized. 

3.7  Flight  Tests 

Data  collection  comprised  of  a  series  of  flight  tests  over  the  pseudolite  array. 
The  flight  test  were  accomplished  over  the  test-range  at  Holloman  AFB  where  the 
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pseudolites  locations  were  established.  The  locations  of  the  pseudolites  based  on  a 
local-level  frame  as  pictured  below  in  Figure  3.4. 


Figure  3.4.  Loclite  Array  Over  Flight  Range 


UHARS  and  CRS  Pallets. 

The  CRS  pallet  is  the  current  legacy  pallet  and  consists  of  a  Honeywell  EEGI, 
Honeywell  SNU  and  Trimble  GPS  receiver.  The  SNU  is  a  navigation  grade  INS, 
and  the  EEGI  achieves  results  higher  than  that  of  a  tactical  grade  inertial.  This 
configuration  is  proven  to  get  decimeter  level  accuracy  for  flight  testing.  The  UHARS 
pallet  consists  of  a  Honeywel  EEGI,  Trimble  GPS  receiver  and  the  Locata  Cooper 
antenna  pseudolite  receiver.  A  series  of  flight  tests  were  accomplished  in  order  to 
ensure  multiple  sets  of  data  were  collected. 

3.8  Summary 

This  section  summarized  the  algorithm  and  application  of  the  measurement  mod¬ 
els  to  obtain  the  pseudolite  navigation  solution.  The  application  of  the  developed 
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algorithms  tropospheric  delay  models  and  cycle  slip  detection  are  a  substantial  im¬ 
portant  to  this  research. 
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IV.  Navigation  Filter  Results  and  Analysis 


4.1  Introduction 

This  chapter  presents  the  results  and  analysis  of  the  navigation  filter  employed 
for  processing  the  data.  The  first  section  of  this  chatper  discusses  carrier-phased 
measurements  from  the  pseudolites  and  the  filters  cycle  slip  detection  and  reset  al¬ 
gorithm.  The  next  section  analyzes  the  effectiveness  of  various  tropospheric  delay 
models  implemented  during  the  course  of  developing  the  filter  navigation  solution. 
The  following  section  discusses  the  navigation  filter  results  from  a  variety  of  tests. 
Finally,  a  summary  section  closes  this  chapter. 

4.2  Carrier-Phase  Measurement  Cycle  Slip  Detection  and  Repair 

Carrier-phase  measurements  provide  the  greatest  level  of  accuracy,  but  they  are 
prone  to  cycle  slips  as  previously  discussed  in  Chapter  3.  Cycle  slips  were  catego¬ 
rized  in  one  of  two  ways  as  previously  discussed  in  Chapter  3:  time  base  cycle  slips 
and  error  magnitude  threshold  cycle  slips.  The  cycle  slips  from  exceeding  the  set 
threshold  accounted  for  the  majority  of  the  slips.  From  Equation  3.20,  the  cycle-slip 
detection  implementation  worked  without  ever  misidentifying  a  false-positive.  Using 
slip  thresholds,  as  defined  in  Equation  3.20  Chapter  3,  of  1.5,  and  3  wavelengths  are 
later  used  for  comparison  in  the  developed  tests.  The  wavelength  of  the  frequency 
used  is  0.1231  meters;  therefore,  the  slip  thresholds  are  equivalent  to  18.47,  and  36.94 
centimeters.  These  thresholds  were  chosen  to  provide  a  constraint  on  the  error  to 
achieve  a  desired  level  of  accuracy,  while  allowing  for  errors  from  other  sources  to 
include  tropospheric  error,  clock  errors,  and  unmodeled  lever  arm  errors.  Using  too 
large  of  a  slip  threshold,  as  in  the  case  of  5  cycles,  or  61.57  centimeters,  lead  to  a 
larger  position  error.  Figures  4.1  and  4.2  show  the  position  error  of  a  5-cycle  and 
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1.5-cycle  threshold  respectivly,  while  keeping  all  other  factors  such  as  lever  arms,  and 
tropospheric  delay  model  used  constant. 


North  Position  Error  (m) 
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Figure  4.1.  Position  Error  when  a  5  Cycle-Slip  Threshold  is  applied 


North  Position  Error  (m) 
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Figure  4.2.  Position  Error  when  a  1.5  Cycle-Slip  Threshold  is  applied 


Figure  4.12  shows  an  overall  better  position  error  when  compared  against  Figure 
4.1  through  visual  inspection.  Further,  for  this  test  a  5-cycle  threhold  compared 
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against  a  3-cycle  threshold  appeared  very  similar  with  outages  occuring  during  the 
same  times  shown  in  Figures  4.3  and  4.4. 

The  change  in  the  slip  threshold,  the  number  of  slips  varied  from  2502  in  the  1.5 
cycle  threshold  case  to  1550  in  the  3  cycle  threshold  case,  during  the  sample  time 
period,  however,  the  slips  are  generally  located  during  the  peroids  where  the  aircraft 
is  banking  as  demonstrated  in  Figures  4.3,  and  4.4. 


PRN  Availability  vs  Time 
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Figure  4.3.  The  LRPN  Availability  Showing  Slip  for  a  1.5  Cycle-Slip  Threshold 
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Figure  4.4.  The  LRPN  Availability  Showing  Slip  for  a  5  Cycle-Slip  Threshold 

This  implies  that  using  larger  slip  thresholds  led  to  higher  initial  position  errors 
and  does  little  to  affect  the  impending  cycle  slips  as  the  receiver  loses  lock  on  the 
pseudolite  signal.  Figures  4.3  and  4.4,  with  a  slip  threshold  set  to  5  cycles  and  1.5 
cycles  respectively,  demonstrates  that  position  error  significantly  decreases  with  a 
higher  cycle  slip  threshold,  as  does  the  uncertainty.  This  indicates  that  the  filter 
is  not  detecting  some  cycle  slips  and  thus  allows  the  introduction  of  error  into  the 
position  measurements. 

The  correlation  between  the  aircrafts  roll  and  cycle  slips  is  most  likely  attributable 
to  multipath  caused  from  signals  reflecting  on  the  aircrafts  fuselage  or  an  error  in  the 
applied  lever  arms.  Figure  4.5,  plots  the  bank  angle  and  the  occurance  of  cycle  slips. 
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Roll  over  Time 


Figure  4.5.  The  Roll  Angle  of  the  Aircraft  with  Cycle  Slips 

Figure  4.5  shows  that  although  cycle  slips  do  not  occur  only  during  a  banking 
maneuver;  however,  they  are  most  prevalent  during  the  peroids  where  the  aircraft 
banks  greater  than  18  degrees.  Further,  since  the  banks  are  at  the  fringe  of  the 
range,  the  dilution  of  precision  is  higher  as  further  discussed  in  Section  4.4. 

The  cycle  slip  repair  process  was  implemented  effectively.  When  examining  the 
carrier  phase  ambiguities,  the  cycle  repair  process  reset  the  state  ambiguity  to  a  new 
value,  as  well  as  resetting  all  cross  correlation  terms  to  zero.  Figure  4.6  plots  the 
change  in  the  ambiguity  estimates  from  the  initial  value  after  the  cycle  slip  reset 
process  against  time.  Cycle  slips  that  exceeded  a  time  threshold  and  cycle-slips  that 
exceeded  a  cycle  threshold  are  shown  along  with  the  ambiguity  estimates  change  from 
the  intial  values. 
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LPRN  100  Ambuguity  Estimates  Change  From  Initialization  vs  Time 


Figure  4.6.  LPRN  100  Cycle  Slips  and  Ambiguity  Tracking 

Figure  4.6  show  the  cycle  slip  repair  process’s  effect  on  the  ambiguities.  The 
change  from  initial  values  of  the  ambiguities  can  change  over  several  cycles;  however, 
when  considering  the  10Hz  data  rate,  each  epoch  change  typically  is  only  a  tenth  of 
a  cycle. 

4.3  Tropospheric  Modeling 

The  tropospheric  delay  models  used  were  a  ray-traced  approach  based  on  the 
modified  Hopheld  model  (RT  mHop),  a  mean  error  correction  based  on  calculating 
the  index  of  refraction  at  the  transmitter  (IRM),  the  Bouska-Raquet  model  (BR), 
and  lastly  a  case  where  no  tropospheric  delay  model  is  applied.  Figure  4.7  plots  the 
tropospheric  delay  correction  that  is  applied  to  LPRN  1400  for  each  of  the  tropospheric 
models  discussed. 
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Troposphere  Error  correction 


Figure  4.7.  Tropospheric  Delay  correction  Applied  to  LRPN  14  over  Sample  Period 

Figure  4.7  shows  that  the  ray-traced  approach  of  the  modified  Hopfield  model 
and  the  Bouska-Raquet  model  bear  similar  tropospheric  delay  model  corrections. 
The  calculated  index  of  refraction  method  for  pseudolites  produce  the  same  general 
shape  to  that  of  the  RT  mHop  approach  and  the  BR  approach;  however,  has  a  flatter 
change  in  magnitude  when  calculating  varying  distances.  From  analysis  in  Chapter 
3,  the  ray-tracing  modified  Hopfield  model  (RT  mHop)  and  the  index  of  refractivity 
model  (IRM)  proved  to  work  the  best  of  the  three  models  when  calculated  using  the 
raw  measurements.  The  Figure  4.7  shows  the  troposphere  corrections  that  were  used 
in  the  sample  period  that  were  applied  to  LPRN  1400. 

4.4  Navigation  Solution  Tests 

Flight  Test  Set-up. 

Data  collected  during  the  flight  test  for  the  UHARS  consisted  of  flying  a  criss¬ 
crossing  pattern  over  the  test  range,  shown  in  Figure  4.8.  The  aircraft  attempted 
to  make  the  turns  not  to  exceed  a  roll  angle  of  25  degrees  to  minimize  the  aircraft 


47 


masking  as  discussed  above  in  Section  4.2.  The  roll  of  the  aircraft  during  the  flight 
over  the  range  is  shown  in  Figure  4.9. 


3D  Flight  Path  in  Lat-Lon  in  Degrees 

- Flight  Trajectory 
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Figure  4.8.  3-D  Flight  Trajectory  Sample  Trajectory  With  Pseudolites 


Figure  4.9.  Aircraft  Roll  Angle  and  Cycle  Slip  vs  Time 


As  previously  discussed,  the  UHARS  EEGI  INS  is  as  the  inertial.  Constant  biases 
are  evident  in  the  1MU  plot  which  is  inductive  of  an  1MU  misalignment  [27]. 


48 


Position  Error  (m) 


Figure  4.10.  Position  Error  of  IMU  Propagation 

From  Figure  4.10,  the  values  for  a  tactical-grade,  IMU  preform  adequately  for 
a  tactical-grade  inertial [2], Evidence  of  the  down  direction  drifting  past  the  accepted 
values  appear  to  be  evident  over  time.  This  is  likely  caused  by  a  misalignment  of  the 
IMU  or  an  unmodeled  lever  arm. 

Developed  Tests  . 

In  order  to  provide  a  comparison  of  the  effectiveness  of  the  developed  algorithm’s 
measurement  model,  cycle  slip  detection  and  repair  process  and  tropospheric  delay 
correction,  six  test  were  developed.  Lever  arms  are  applied  from  the  EEGI  on  the 
UHARS  rack  to  the  Cooper  antenna,  which  receive  dthe  Locata  pseudolites  signals. 
The  Liters  initial  position,  velocity,  and  attitude  is  initialized  from  the  truth  reference 
for  that  instant  in  time.  The  truth  reference  that  is  used  to  compare  the  filter  solution 
was  the  746Th  Test  Squadron  blended  solution  that  also  originates  at  the  EEGI  on  the 
UHARS  rack.  This  allows  the  Liter  generated  solution  and  truth  reference  solution 
to  have  the  same  origin  and  therefore  no  additional  lever  arms  are  required  to  be 
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applied.  The  tropospheric  delay  models,  slip  thresholds,  and  number  of  pseudorange 
updates  recorded  after  a  cycle  slip,  as  described  in  Section  3.6  in  Chapter  III,  vary 
in  each  test.  A  breakdown  of  the  tests  are  in  Table  4.1  with  varying  tropospheric 
models,  cycle  slip  thresholds,  and  pseudorange  updates  after  a  cycle  slip  occurs. 

Table  4.1.  Developed  Tests  for  Comparison 


Scenarios 


Scenario  1 

Trop.  Model 

None 

cycles 

3 

PR  updates 

2 

Scenario  2 

RT  mHop 

1.5 

2 

Scenario  3 

RT  mHop 

3 

2 

Scenario  4 

RT  mHop 

3 

0 

Scenario  5 

RT  mHop 

3 

oo 

Scenario  6 

IRM 

3 

2 

The  term  RT  mHop  refers  to  the  ray-tracing  method  of  the  modified  Hopheld 
model,  and  the  term  IRA4,  refers  to  the  index  of  refraetivity  model  developed.  Since 
the  tropospheric  correction  between  the  RT  mHop  and  the  BR  model  were  close, 
only  the  RT  mHop  model  was  used  in  the  tests  developed.  The  sample  period  chosen 
for  the  flight  tests  were  North-South  flights  shown  in  Figure  4.11.  The  blue  circles 
represent  the  pseudolite  locations,  and  the  red  line  shows  the  flight  path  that  begins 
on  the  bottom  right  of  the  plot  and  continues  toward  the  left.  For  each  test,  the 
position  error  and  measurement  residuals  are  plotted  for  comparison. 


50 


Figure  4.11.  2-D  Flight  Trajectory  During  Test  Peroid  with  Pseudolites 

Analysis  of  Developed  Tests. 

After  the  six  tests  are  established  above,  the  position  results  are  plotted  in  this  sec¬ 
tion,  analyzed  and  compared.  Further,  the  filters  residuals  are  plotted  and  compared. 
Together,  plots  proved  a  measure  of  the  navigation  filter  effectiveness. 

Test  1. 

The  first  scenario,  does  not  apply  a  tropospheric  model,  and  uses  two  pseudorange 
updates  after  a  cycle-slip  with  a  cycle-slip  threshold  set  to  3  cycles.  Figure  4.12,  plots 
Test  l’s  position  error. 
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Figure  4.12.  Position  Error  of  Test  1.  No  tropospheric  delay  model  applied,  with  a  slip 
threshold  of  3  cycles  and  2  pseudorange  updates  upon  clcyle  slip  reset. 

When  examining  the  position  error  upon  initialization,  in  all  the  North,  East  and 
down  directions  the  error  rapidly  decreases  as  the  filter  converges  to  a  navigation 
solution.  In  the  north  direction  that  the  error  continues  to  exceed  the  uncertainty 
expected;  however,  the  error’s  magnitude  decreases  over  time,  and  converges  to  an 
error  value  within  the  sigma  values.  The  east  position  error  begins  to  drift  as  time 
continues.  Lastly,  the  down  direction  error  appears  constant,  only  slightly  improving 
with  time,  appearing  to  have  a  near  constant  bias.  The  error  is  most  likely  attributable 
to  tropospheric  error  model  and/or  unmodeled  lever  arms.  Some  of  the  tropospheric 
error  in  this  test  are  absorbed  into  other  bias  error  terms.  Figure  4.13  shows  the 
carrier  phase  residuals  of  test  1. 
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Figure  4.13.  Carrier  Phase  Measurement  Residuals  of  Test  1 

Test  2. 

The  filter  residual  errors  peaked  during  the  cycle  clip  repair  process.  Somewhat 
higher  residuals  were  expected,  based  on  the  position  error,  which  indicated  that  the 
filter  converged  onto  a  solution  with  inherent  error. 

In  test  2,  the  modified  Hopfield  tropospheric  delay  model  was  added  with  again 
two  pseudorange  measurements  added  after  a  cycle  slip,  and  a  cycle  slip  threshold  of 
1.5  cycles  is  set.  The  position  error  is  shown  in  4.14. 
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Figure  4.14.  Position  Error  of  Test  2.  RT  mHoptropospheric  delay  model  applied,  with 
a  slip  threshold  of  1.5  cycles  and  2  pseudorange  updates  upon  clcyle  slip  reset. 

The  position  error  general  magnitudes  in  all  directions  are  relatively  consistent. 
Some  indication  of  an  additional  lever  arm  that  is  not  modeled  appears  to  exist  as 
the  error  in  the  east  direction  is  affected  by  the  direction  of  flight.  The  carrier  phase 
residuals  are  plotted  in  4.15. 


Resids 


Figure  4.15.  Carrier  Phase  Measurement  Residuals  of  Test  2 
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The  filter  residuals  look  very  similar  to  that  of  the  first  scenario,  despite  a  tropo¬ 
spheric  delay  model  being  added.  Figure  4.15  show  that  the  residuals  spike  during 
cycle  slips,  and  quickly  return  to  low  (less  than  2)  values. 

Test  3. 

Test  3  is  much  like  Test  2,  however,  a  looser  slip  threshold  of  3  cycles  is  applied. 
The  position  error  is  shown  in  Figure  4.16. 
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Figure  4.16.  Position  Error  of  Test  3.  RT  mHop  tropospheric  delay  model  applied, 
with  a  slip  threshold  of  3  cycles  and  2  pseudorange  updates  upon  clcyle  slip  reset. 

The  looser  slip  threshold  overall  appear  to  have  small  effect  on  the  position  error. 
Without  further  statistical  analysis,  little  in  inferred  about  the  position  error.  Figure 
4.17  shows  the  carrier  phase  measurement  residuals. 
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Figure  4.17.  Carrier  Phase  Measurement  Residuals  of  Test  3 

Upon  comparison  of  the  carrier  phase  residuals  of  Figure  4.15  and  Figure  4.17 
insignificant  change  is  visible  seen.  This  infers  that  the  change  of  a  slip  threshold 
from  1.5  cycles  to  3  cycles  has  a  small  effect  on  the  overall  solution  accuracy  and 
performance. 

Test  4. 

Test  4  also  closely  follows  Test  2,  with  the  exception  that  this  test  does  not  incor¬ 
porate  any  pseudorange  measurements.  The  strictly  carrier  phase  position  solution 
is  shown  in  Figure  4.18. 
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Figure  4.18.  Position  Error  of  Test  4.  RT  mHop  tropospheric  delay  model  applied, 
with  a  slip  threshold  of  3  cycles  and  no  pseudorange  updates  upon  clcyle  slip  reset. 

The  position  solution  for  this  test  had  smaller  sigma  values  since  the  carrier-phase 
sigma  value  was  set  to  0.07  meters  and  a  slightly  worst  position  estimate  as  shown 
in  Figure  4.18  when  compared  against  Figure  4.14.  The  measurement  residuals  are 
shown  in  Figure  4.19  for  this  test. 


Phase  Residuals 


Figure  4.19.  Carrier  Phase  Measurement  Residuals  of  Test  4 
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Figure  4.19  again  to  have  the  larger  residual  errors  during  cycle  slips  then  converge 
to  low  values  (less  than  2)  as  in  previous  test  cases. 

Test  5. 

Test  5,  much  like  Test  4,  closely  follows  Test  2,  however,  incorporates  every  avail¬ 
able  pseudorange  measurement.  The  postion  error  of  Test  5  is  shown  in  Figure  4.20. 


tgps  (sec) 


Figure  4.20.  Position  Error  of  Test  5.  RT  mHop  tropospheric  delay  model  applied, 
with  a  slip  threshold  of  3  cycles  and  an  infinite  number  of  pseudorange  updates  upon 
clcyle  slip  reset. 

Since  the  pseudorange  measurements  are  noisy,  a  higher  uncertainty  (4  meters) 
is  assigned  to  them.  The  filters  position  solution  stayed  within  the  higher  bounds, 
although  the  overall  accuracy  was  not  improved.  Figure  4.21  shows  the  pseudorange 
residuals  and  Figure  4.22  shows  only  the  carrier-phase  residuals. 
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Figure  4.21.  Pseudorange  Residuals  of  Test  5.  IRM  tropospheric  delay  model  applied, 
with  a  slip  threshold  of  3  cycles  and  unlimited  pseudorange  updates  upon  clcyle  slip 
reset. 


Phase  Residuals 


Figure  4.22.  Phase  Residuals  of  Test  5.  IRM  tropospheric  delay  model  applied,  with  a 
slip  threshold  of  3  cycles  and  unlimited  pseudorange  updates  upon  clcyle  slip  reset. 

The  carrier-phase  residuals  shown  in  Figure  4.22  follow  very  closely  to  the  previous 
test  residuals.  The  pseudorange  residuals  are  much  higher  than  the  carrier-phase 
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residuals  as  shown  in  4.21  due  to  the  noisy  nature  of  pseudorange  measurements. 


Test  6. 

Lastly,  Test  6  closely  followed  Test  2  using  a  cycle  threshold  of  3  cycles,  and 
incorporating  two  pseudoranges  after  a  cycle  slip;  however,  this  model  uses  the  index 
of  refrativity  tropospheric  delay  model.  .  Figure  4.23  plots  the  position  error  for  this 
test. 
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Figure  4.23.  Position  Error  of  Test  6.  IRM  tropospheric  delay  model  applied,  with  a 
slip  threshold  of  3  cycles  and  two  pseudorange  updates  upon  clcyle  slip  reset. 

Test’s  6  appeared  to  have  the  much  worst  overall  results  in  terms  of  position 
error  visibly  seen  in  Figure  4.23.  This  indicates  that  the  tropospheric  delay  model 
preformed  less  well  than  that  of  the  ray-tracing  modified  Hopheld  model.  The  carrier- 
phase  measurement  residuals  shown  in  Figure  4.24. 
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Res  ids 


Figure  4.24.  Carrier  Phase  Measurement  Residuals  of  Test  6 

Although  the  position  error  was  significantly  higher  as  the  previous  tests,  the 
carrier  phase  measurement  residuals  appeared  only  slightly  worst  (higher). 

Summary  of  Tests. 

From  analysis,  a  level  arm  appears  to  be  unaccounted  for  between  the  navigation 
solution  and  the  truth.  This  is  evident  from  the  correlation  between  the  direction 
of  error  in  the  position  error  plot  and  the  heading  of  the  aircraft.  Error  that  exist 
from  a  lever  arm  does  not  account  for  all  the  error  that  exists.  As  the  aircraft  travels 
along  the  outer  edge  of  the  Loclite  array,  the  error  grows  in  magnitude.  This  trend 
is  visible  from  all  six  tests  performed. 
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2DRMS  and  3DRMS  Analysis. 


The  2D  and  3DRMS  are  good  metrics  for  postion  accuracy [15].  The  2DRMS  uses 
the  north  and  east  position  errors  only  are  are  calculated  by  using  4.1.  the  3DRMS 
incorporates  a  3-axis  postion  accuracy  metric  and  is  calculate  using  Equation  4.2. 
The  5Pn,  the  5Pn,  and  the  8Pn  terms  represent  the  position  errors  in  north,  east  and 
down  direction  respectively. 


2D  -  RMS  =  y/6P%  +  8PI  (4.1) 

3 D  -  RMS  =  yJsP*  +  6 Pi  +  SPt |  (4.2) 

The  mean  of  the  tests  were  were  also  compared  againt  a  single  pass  in  the  center  of 
the  Loclite  array.  This  single  pass  took  place  over  a  200  second  period  and  is  plotted 
below  in  Figure  4.25. 
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Figure  4.25.  Center  Pass  Flight  Trajectory 


The  2D  and  3DRMS  were  compared  to  one  another  taking  the  mean  of  each  test 
and  during  the  center  pass  over  a  200  second  period. 
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Table  4.2.  2D  and  3D  Root  Mean  Square  Errors 


2DRMS 

3D  RMS 

Test 

Center  Pass 

(m)  Period  Mean  (m) 

Center  Pass  (m) 

Period  Mean  (m 

1 

0.50 

0.74 

0.74 

1.63 

2 

0.26 

0.74 

0.59 

2.05 

3 

0.21 

0.59 

0.57 

1.51 

4 

0.23 

0.2852 

0.62 

0.55 

5 

0.24 

.0.61 

0.60 

1.61 

6 

0.31 

0.71 

0.71 

1.97 

7 

0.18 

0.60 

0.45 

1.36 

Dilution  of  Precision  -  DOP. 


The  dilution  of  precision  (DOP)  is  analyzed  to  estimate  how  much  the  measure¬ 
ment  error  will  affect  the  final  state  estimation.  The  DOP  is  simply  a  magnifying 
effect  based  on  the  geometry  of  the  transmitters.  More  accurately,  the  DOP  is  based 
on  the  covariance  matrices  and  the  clock  error.  To  find  the  DOP  the  user  to  pseudo- 
lite  the  equation  is  first  linearized  to  a  nominal  user  position  as  defined  in  Equation 


4.3. 


rp  — 
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zu0 


p 

PuO 


Azu  +  cAT 


(4.3) 


In  the  above  equation,  the  p  term  represents  the  pseudolite,  where  the  u  term 
refers  to  the  user  in  the  local  level  frame,  the  0  subscript  represents  the  nominal 
position,  r  is  the  range  respect  to  the  nominal  user  position,  x ,  y,  z  are  axis  positions 
and  cAtu  is  the  clock  error  in  units  of  meters.  The  linearized  ranging  equation  is 
shown  in  equation  4.3  and  is  rearranged  to  have  the  form  A p  =  HAx„.  The  DOP  is 
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preformed  by  calculating  the  variances  by  preforming  D  =  (H2H)_1.  The  diagonals 
represents  the  covariance  of  one  axis  D n  =  cby  D2 2  =  cr^,  and  D33  =  cr^. 

The  variances  required  for  calculating  the  DOP  is  obtained  from  the  this  matrix 
D  =  (H2H)_1.  Each  diagonal  represents  the  covariance  of  one  axis  Du  =  a 
D22  =  CTy,  etc.  The  equation  for  the  different  types  of  DOP  that  were  implemented 
in  this  research  are  shown  in  4.4. 

Position  DOP 
=  \J  Du  +  D22  +  D33 

DOP  =  (4.4) 

Horizontal  DOP 

=  y/  Du  +  D22 

The  DOP  values  were  calculate  at  the  center  of  the  range,  the  northern  fringe  of 
the  range,  and  over  one  of  the  tests.  Figure  4.26  depicts  where  the  DOP  measurements 
at  the  Northern  fringe  and  center  postions. 


Loclite  Array 


Easting  kilometers 


Figure  4.26.  Flight  trajectory  showing  positions  where  DOP  was  measured 


Table  4.3  show  the  sample  DOP  values  taken  at  the  center,  norther  edge  of  the 
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array,  and  the  period  mean.  Figure  4.26  shows  the  DOP  values  over  the  test  period. 


Table  4.3.  Sample  Dilution  of  Precision  Values 


Calculated  Dilution  of  Precision 


Center 

N 

E 

S 

W 

Mean 

HDOP 

0.62 

1.89 

1.67 

1.721 

1.62 

1.43 

PDOP 

1.42 

3.13 

2.97 

2.95 

2.98 

2.58 

4.5  Summary 

This  chapter  discussed  the  navigation  filter  parameters,  description  of  the  test 
developed  and  result.  This  research  is  the  first  time  that  Locata  pseudolites  were 
tightly-coupled  with  an  1MU  in  an  aircraft  scenario.  No  single  test  proved  to  obtain 
the  desired  results  of  decimeter  level  position  accuracy  which  was  indicative  of  un¬ 
modeled  errors.  This  leads  to  the  conclusion  that  unmodeled  errors  exist  preventing 
the  desired  decimeter  level  of  accuracy  to  be  achieved. 
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V.  Conclusions  and  Recommendations 


5.1  Conclusions 

This  thesis  presented  the  initial  implementation  of  a  tightly-coupled  Locata  pseu- 
dolite/INS  for  developing  a  navigation  solution  for  the  UHARS  program.  The  objec¬ 
tive  of  decimeter-level  accuracy  navigation  solution  results  was  not  obtained  for  any 
the  test  cases  performed.  Sub-meter  accuracy,  however,  was  consistently  achieved 
which  indicates  unmodeled  errors  still  exist  within  the  solution.  The  tropospheric 
delay  model,  unmodeled  lever  arms,  and  filter  tuning  are  likely  causes  of  the  solution 
not  to  receive  the  desired  level  of  accuracy  expected  for  a  carrier-phase  measurements 
solution. 

A  cycle  slip  detection  and  repair  process  was  implemented  for  detecting  failures 
in  pseudolite  phase  measurements.  The  pseudolite  suffered  from  a  high  rate  of  cycle 
slips  and  the  implemented  cycle  slip  detection  and  repair  process  appeared  to  work 
well.  The  ambiguities  states  under  the  tests  parameters  provide  tracked  well  give  the 
parameters.  Using  a  cycle  slip  threshold  between  1.5  to  3  cycles  both  proved  to  lead 
to  submeter-level  accuracy.  In  Chapter  IV  when  the  slip  threshold  was  extended  to 
5  cycles,  the  position  error  increased  as  error  was  introduced  into  the  filter  by  the 
looser  threshold  value. 

A  further  exploration  of  a  tropospheric  models  can  decrease  the  position  error  as 
well.  The  index  of  refraction  method  faired  worse  than  that  of  the  ray-traced  modified 
Hopfield  model  and  the  Bouska-Raquest  model  whose  results  closely  followed  that  of 
the  the  ray-traced  modified  Hopfield  model. 

Since  there  appeared  to  be  constant  biases  in  the  position  error  depending  on  the 
flight  heading,  unmodeled  lever  arms  were  the  likely  cause.  Further,  this  appeared  as 
a  reoccuring  source  of  error  for  every  test  case  developed,  thus  making  the  possibility 
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unlikely  that  other  error  factors  were  leading  to  the  constant  biases. 


Lastly,  the  carrier-phase  solution  with  only  a  few  pseudorange  updates  preformed 
better  than  that  of  the  test  case  where  every  pseudorange  available  was  used. 

5.2  Recommendations  for  Future  Research 

There  exist  several  recommendations  for  future  research.  Each  recommendation 
will  further  the  advance  the  research  and  quality  of  the  navigation  solution 

•  The  ray-tracing  modified  Hopheld  tropospheric  delay  model,  as  described  in 
Chapter  III,  may  have  preformed  better  if  data  was  available  for  temperature, 
pressure,  and  relative  humidity  changes  were  recorded  over  the  path  of  the  signal 
transmission,  rather  than  only  available  at  the  transmitters  location.  Given  the 
change  in  temperature,  pressure,  and  relativity,  a  more  accurate  tropospheric 
delay  model  can  be  achieved. 

•  Implementing  software  for  the  Kalman  filter  to  accept  dual  frequencies.  The 
implemented  filter  in  this  research  only  used  one  of  the  two  transmitted  fre¬ 
quencies.  By  using  both  frequencies,  a  more  robust  cycle  slip  algorithm  can  be 
developed.  Further,  using  dual  frequencies  can  lead  to  several  other  opportuni¬ 
ties  to  increase  filter  accuracy. 

•  Filtertuning  should  be  continued  in  order  to  improve  filter  performance. 

•  Lever  arms  should  be  either  re-measured,  or,  alternatively,  lever  arms  could  be 
estimated  in  the  filter  in  order  to  reduce  errors  due  to  lever  arms  error. 
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